home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Language/OS - Multiplatform Resource Library
/
LANGUAGE OS.iso
/
cpp_libs
/
rwvector.lha
/
RWVector2.1
/
src
/
dcosine.cc
< prev
next >
Wrap
C/C++ Source or Header
|
1989-08-18
|
4KB
|
168 lines
/*
* Definitions for the double precision cosine server
*
* Copyright (C) 1988, 1989.
*
* Dr. Thomas Keffer
* Rogue Wave Associates
* P.O. Box 85341
* Seattle WA 98145-1341
*
* Permission to use, copy, modify, and distribute this
* software and its documentation for any purpose and
* without fee is hereby granted, provided that the
* above copyright notice appear in all copies and that
* both that copyright notice and this permission notice
* appear in supporting documentation.
*
* This software is provided "as is" without any
* expressed or implied warranty.
*
*
* @(#)dcosine.cc 2.1 8/18/89
*/
#include "rw/DCosineFFT.h"
#include "rw/mathpack.h"
DoubleCosineServer::DoubleCosineServer()
{
server_N=0;
}
DoubleCosineServer::DoubleCosineServer(unsigned N)
{
setOrder(N);
}
DoubleCosineServer::DoubleCosineServer(const DoubleCosineServer& s)
: (s), sins(s.sins)
{
server_N = s.server_N;
}
void
DoubleCosineServer::operator=(const DoubleCosineServer& s)
{
DoubleFFTServer::operator=(s);
server_N = s.server_N;
sins.reference(s.sins);
}
void
DoubleCosineServer::setOrder(unsigned N)
{
server_N = N;
DoubleFFTServer::setOrder(N/2+1);
sins.reference(2.0*sin(DoubleVec(N-1,M_PI/N, M_PI/N)));
}
/*
The fourier transform of a real even sequence is a cosine transform.
The result is also a real even sequence. This is prodedure #6 from
Cooley, et al. (1970) J. Sound Vib. (12) 315--337.
*/
DoubleVec
DoubleCosineServer::cosine(const DoubleVec& v)
{
unsigned Nstore = v.length();
unsigned N = Nstore-1;
unsigned Nm1 = N-1;
unsigned Nhalf = N/2;
unsigned Nhalfp1 = Nhalf+1;
unsigned Nhalfm1 = Nhalf-1;
// Various things that could go wrong:
checkOdd(Nstore);
if( N != server_N )setOrder(N);
/* The results will appear here, and will consist of length() points.*/
DoubleVec results(Nstore);
/* X will be a complex conjugate even vector */
DComplexVec X(Nhalfp1);
X.slice(1,Nhalfm1,1) =
DComplexVec(v.slice(2,Nhalfm1,2), // Real part
v.slice(1,Nhalfm1,2) - v.slice(3,Nhalfm1,2)); // Imag part
X(0) = DComplex(v(0));
X(Nhalf) = DComplex(v(N));
/* Take the IDFT. The transform of a complex conjugate even vector
N/2+1 points long is a real vector, N points long. */
DoubleVec A = DoubleFFTServer::ifourier(X);
DoubleVec temp1 = A.slice(1,Nm1,1);
DoubleVec temp2 = A.slice(Nm1,Nm1,-1);
results.slice(1,Nm1,1) = 0.5*( (temp1 + temp2) - (temp1 - temp2)/sins);
/* Special procedure required for the first and last points. */
// Sum all of the odd points. Count all but v(N) twice.
double A2 = 2.0*sum(v.slice(1,Nhalf,2));
if( N%2) A2 += v(N); // If N is odd, include v(N)
results(0) = A(0) + A2;
results(N) = A(0) - A2;
return results;
}
/*
The fourier transform of a real odd sequence is a sine transform.
The result is also a real odd sequence. This is prodedure #7 from
Cooley, et al. (1970) J. Sound Vib. (12) 315--337.
*/
DoubleVec
DoubleCosineServer::sine(const DoubleVec& v)
{
unsigned Nstore = v.length();
unsigned N = Nstore+1;
unsigned Nm1 = N-1;
unsigned Nhalf = N/2;
unsigned Nhalfp1 = Nhalf+1;
unsigned Nhalfm1 = Nhalf-1;
// Various things that could go wrong:
checkOdd(Nstore);
if(N != server_N)setOrder(N);
/* X will be a complex conjugate even vector */
DComplexVec X(Nhalfp1);
X.slice(1,Nhalfm1,1) =
DComplexVec(v.slice(0,Nhalfm1,2) - v.slice(2,Nhalfm1,2), // real part
-v.slice(1,Nhalfm1,2)); // imag part
X(0) = DComplex(-2.0*v(0));
X(Nhalf) = DComplex( 2.0*v(N-2));
/* Take the IDFT. The transform of a complex conjugate even vector
N/2+1 points long is a real vector, N points long. */
DoubleVec A = DoubleFFTServer::ifourier(X);
DoubleVec temp1 = A.slice(1,Nm1,1);
DoubleVec temp2 = A.slice(Nm1,Nm1,-1);
DoubleVec results = 0.5*( (temp1 - temp2) - (temp1 + temp2)/sins);
return results;
}
void
DoubleCosineServer::checkOdd(int l)
{
if( !(l%2) ){
char msg[80];
sprintf(msg, "Length (%d) of a real series must be odd.", l);
RWnote("DoubleCosineServer::[i]fourier()", msg);
RWerror(DEFAULT);
}
}